/*
 * Copyright 2010-2012 Susanta Tewari. <freecode4susant@users.sourceforge.net>
 *
 * This program is free software: you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program.  If not, see <http://www.gnu.org/licenses/>.
 */

package bd.org.apache.commons.math.util;

import bd.org.apache.commons.math.exception.ConvergenceException;
import bd.org.apache.commons.math.exception.MaxCountExceededException;
import bd.org.apache.commons.math.exception.util.LocalizedFormats;

/**
 * Provides a generic means to evaluate continued fractions.  Subclasses simply
 * provided the a and b coefficients to evaluate the continued fraction.
 * <p>
 * References:
 * <ul>
 * <li><a href="http://mathworld.wolfram.com/ContinuedFraction.html">
 * Continued Fraction</a></li>
 * </ul>
 * </p>
 *
 * @version $Id: ContinuedFraction.java 1244107 2012-02-14 16:17:55Z erans $
 */
public abstract class ContinuedFraction {

    /**
     * Maximum allowed numerical error.
     */
    private static final double DEFAULT_EPSILON = 10e-9;

    /**
     * Default constructor.
     */
    protected ContinuedFraction() {

        super();
    }

    /**
     * Access the n-th a coefficient of the continued fraction.  Since a can be
     * a function of the evaluation point, x, that is passed in as well.
     *
     * @param n the coefficient index to retrieve.
     * @param x the evaluation point.
     * @return the n-th a coefficient.
     */
    protected abstract double getA(int n, double x);


    /**
     * Access the n-th b coefficient of the continued fraction.  Since b can be
     * a function of the evaluation point, x, that is passed in as well.
     *
     * @param n the coefficient index to retrieve.
     * @param x the evaluation point.
     * @return the n-th b coefficient.
     */
    protected abstract double getB(int n, double x);


    /**
     * Evaluates the continued fraction at the value x.
     *
     * @param x the evaluation point.
     * @return the value of the continued fraction evaluated at x.
     */
    public double evaluate(double x) {

        return evaluate(x, DEFAULT_EPSILON, Integer.MAX_VALUE);
    }


    /**
     * Evaluates the continued fraction at the value x.
     *
     * @param x the evaluation point.
     * @param epsilon maximum error allowed.
     * @return the value of the continued fraction evaluated at x.
     */
    public double evaluate(double x, double epsilon) {

        return evaluate(x, epsilon, Integer.MAX_VALUE);
    }


    /**
     * Evaluates the continued fraction at the value x.
     *
     * @param x the evaluation point.
     * @param maxIterations maximum number of convergents
     * @return the value of the continued fraction evaluated at x.
     */
    public double evaluate(double x, int maxIterations) {

        return evaluate(x, DEFAULT_EPSILON, maxIterations);
    }


    /**
     * <p>
     * Evaluates the continued fraction at the value x.
     * </p>
     * <p>
     * The implementation of this method is based on equations 14-17 of:
     * <ul>
     * <li>
     * Eric W. Weisstein. "Continued Fraction." From MathWorld--A Wolfram Web
     * Resource. <a target="_blank"
     * href="http://mathworld.wolfram.com/ContinuedFraction.html">
     * http://mathworld.wolfram.com/ContinuedFraction.html</a>
     * </li>
     * </ul>
     * The recurrence relationship defined in those equations can result in
     * very large intermediate results which can result in numerical overflow.
     * As a means to combat these overflow conditions, the intermediate results
     * are scaled whenever they threaten to become numerically unstable.</p>
     *
     * @param x the evaluation point.
     * @param epsilon maximum error allowed.
     * @param maxIterations maximum number of convergents
     * @return the value of the continued fraction evaluated at x.
     */
    public double evaluate(double x, double epsilon, int maxIterations) {

        double p0            = 1.0;
        double p1            = getA(0, x);
        double q0            = 0.0;
        double q1            = 1.0;
        double c             = p1 / q1;
        int    n             = 0;
        double relativeError = Double.MAX_VALUE;

        while ((n < maxIterations) && (relativeError > epsilon)) {

            ++n;

            double  a        = getA(n, x);
            double  b        = getB(n, x);
            double  p2       = a * p1 + b * p0;
            double  q2       = a * q1 + b * q0;
            boolean infinite = false;

            if (Double.isInfinite(p2) || Double.isInfinite(q2)) {

                /*
                 * Need to scale. Try successive powers of the larger of a or b
                 * up to 5th power. Throw ConvergenceException if one or both
                 * of p2, q2 still overflow.
                 */
                double       scaleFactor     = 1d;
                double       lastScaleFactor = 1d;
                final int    maxPower        = 5;
                final double scale           = FastMath.max(a, b);

                if (scale <= 0) {    // Can't scale

                    throw new ConvergenceException(
                        LocalizedFormats.CONTINUED_FRACTION_INFINITY_DIVERGENCE, x);
                }

                infinite = true;

                for (int i = 0; i < maxPower; i++) {

                    lastScaleFactor = scaleFactor;
                    scaleFactor     *= scale;

                    if ((a != 0.0) && (a > b)) {

                        p2 = p1 / lastScaleFactor + (b / scaleFactor * p0);
                        q2 = q1 / lastScaleFactor + (b / scaleFactor * q0);

                    } else if (b != 0) {

                        p2 = (a / scaleFactor * p1) + p0 / lastScaleFactor;
                        q2 = (a / scaleFactor * q1) + q0 / lastScaleFactor;
                    }

                    infinite = Double.isInfinite(p2) || Double.isInfinite(q2);

                    if (!infinite) {
                        break;
                    }
                }
            }

            if (infinite) {

                // Scaling failed
                throw new ConvergenceException(
                    LocalizedFormats.CONTINUED_FRACTION_INFINITY_DIVERGENCE, x);
            }

            double r = p2 / q2;

            if (Double.isNaN(r)) {

                throw new ConvergenceException(LocalizedFormats.CONTINUED_FRACTION_NAN_DIVERGENCE,
                                               x);
            }

            relativeError = FastMath.abs(r / c - 1.0);

            // prepare for next iteration
            c  = p2 / q2;
            p0 = p1;
            p1 = p2;
            q0 = q1;
            q1 = q2;
        }

        if (n >= maxIterations) {

            throw new MaxCountExceededException(LocalizedFormats.NON_CONVERGENT_CONTINUED_FRACTION,
                    maxIterations, x);
        }

        return c;
    }
}
